Assignment 4: System Planning¶
In this group work, we will investigate scenarios for a 100% renewable electricity system for the country ofSlovakia, which is a landlocked, mountainous country in Europe, lies roughly between latitudes 47° and 50° N, and longitudes 16° and 23° E. Forests cover 41% of Slovak land surface. [wikipedia].¶
In [1]:
# importing the necessary libraries & Packages:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
import pypsa
import cartopy.feature as cfeature
import cartopy.crs as ccrs
import cartopy
import warnings
import atlite
import rasterio
from rasterio.plot import show
import xarray as xr
import netCDF4
warnings.filterwarnings('ignore')
from atlite.gis import ExclusionContainer
from atlite.gis import shape_availability
import rasterio as rio
import os
import country_converter as coco
import logging
logging.basicConfig(level=logging.INFO)
from shapely.geometry import Polygon, MultiPolygon, Point
In [2]:
# extracting GADM data belonging to Slovakia & creating the shape of the country including its regions:
svk_regions = gpd.read_file('gadm_410-levels-ADM_1-SVK.gpkg')
svk_regions = svk_regions.set_index('NAME_1')
display(svk_regions)
shape_svk = svk_regions.to_crs(3035).geometry # Slovokia is in Europe (crs = 3035)
shape_svk.plot()
| GID_0 | GID_1 | COUNTRY | geometry | |
|---|---|---|---|---|
| NAME_1 | ||||
| Banskobystrický | SVK | SVK.1_1 | Slovakia | MULTIPOLYGON (((19.45671 48.09119, 19.45659 48... |
| Bratislavský | SVK | SVK.2_1 | Slovakia | MULTIPOLYGON (((17.31282 48.45557, 17.31436 48... |
| Košický | SVK | SVK.3_1 | Slovakia | MULTIPOLYGON (((22.10827 48.38136, 22.10702 48... |
| Nitriansky | SVK | SVK.4_1 | Slovakia | MULTIPOLYGON (((18.83296 47.81751, 18.81611 47... |
| Prešovský | SVK | SVK.5_1 | Slovakia | MULTIPOLYGON (((21.54594 48.79453, 21.54533 48... |
| Trenčiansky | SVK | SVK.6_1 | Slovakia | MULTIPOLYGON (((17.41709 48.83292, 17.42872 48... |
| Trnavský | SVK | SVK.7_1 | Slovakia | MULTIPOLYGON (((17.45999 48.14781, 17.46112 48... |
| Žilinský | SVK | SVK.8_1 | Slovakia | MULTIPOLYGON (((19.02072 48.87508, 19.01839 48... |
Out[2]:
<Axes: >
In [3]:
###### !!! this section is performed only for presentation purposes !! #########
df = pd.read_csv('global_power_plant_database.csv') # Reading the data frame included all power plant data of all countries
df_svk = df[df.iloc[:, 0] == 'SVK'] # extracting the Slovakia-related data
display(df_svk.head(5))
df_total_generation_groupby = df_svk.groupby('primary_fuel').estimated_generation_gwh_2017.sum()
#df_total_generation_groupby.head(5)
# Plotting the pie chart of total generation per technology in Slovakia:
plt.figure(figsize=(8, 5))
plt.pie(df_total_generation_groupby, labels=df_total_generation_groupby.index, autopct='%1.1f%%', startangle=90)
plt.legend(df_total_generation_groupby.index, bbox_to_anchor=(1, 0.5), loc="center left", title="Energy Carriers")
plt.title('Estimated Generation GWh (2017) by Primary Fuel for Slovakia')
plt.savefig('Estimated Generation(2017).png')
plt.show()
print()
print()
geometry = gpd.points_from_xy(df_svk["longitude"], df_svk["latitude"])
gdf_svk = gpd.GeoDataFrame(df_svk, geometry=geometry, crs=4326)
| country | country_long | name | gppd_idnr | capacity_mw | latitude | longitude | primary_fuel | other_fuel1 | other_fuel2 | ... | estimated_generation_gwh_2013 | estimated_generation_gwh_2014 | estimated_generation_gwh_2015 | estimated_generation_gwh_2016 | estimated_generation_gwh_2017 | estimated_generation_note_2013 | estimated_generation_note_2014 | estimated_generation_note_2015 | estimated_generation_note_2016 | estimated_generation_note_2017 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 19861 | SVK | Slovakia | Badín | WKS0063339 | 6.0 | 48.2650 | 19.0560 | Solar | NaN | NaN | ... | 8.29 | 8.91 | 8.90 | 8.58 | 8.36 | SOLAR-V1-NO-AGE | SOLAR-V1-NO-AGE | SOLAR-V1-NO-AGE | SOLAR-V1-NO-AGE | SOLAR-V1-NO-AGE |
| 19862 | SVK | Slovakia | Bohunice Nuclear Power Plant Slovakia | GEODB0002385 | 880.0 | 48.4960 | 17.6816 | Nuclear | NaN | NaN | ... | NaN | NaN | NaN | NaN | 6805.78 | NO-ESTIMATION | NO-ESTIMATION | NO-ESTIMATION | NO-ESTIMATION | CAPACITY-FACTOR-V1 |
| 19863 | SVK | Slovakia | Cierny Vah Pumped Storage Hydroelectric Power ... | GEODB0042635 | 735.0 | 49.0088 | 19.9122 | Hydro | NaN | NaN | ... | 2597.78 | 2210.18 | 2927.73 | 2216.24 | 2210.18 | HYDRO-V1 | HYDRO-V1 | HYDRO-V1 | HYDRO-V1 | HYDRO-V1 |
| 19864 | SVK | Slovakia | Cunovo Hydroelectric Power Plant Slovakia | GEODB0042641 | 24.0 | 48.0299 | 17.2254 | Hydro | NaN | NaN | ... | 118.03 | 117.69 | 129.12 | 122.91 | 130.73 | HYDRO-V1 | HYDRO-V1 | HYDRO-V1 | HYDRO-V1 | HYDRO-V1 |
| 19865 | SVK | Slovakia | Dubnica nad Vahom Hydroelectric Power Plant Sl... | GEODB0042662 | 16.5 | 48.9638 | 18.1454 | Hydro | NaN | NaN | ... | 73.11 | 55.75 | 63.33 | 52.24 | 55.80 | HYDRO-V1 | HYDRO-V1 | HYDRO-V1 | HYDRO-V1 | HYDRO-V1 |
5 rows × 36 columns
Out[3]:
<GeoAxes: >
In [4]:
# Reading the file included all of the Exclusive economic zones of the countries, which is not the case for Slovakia because it is a landlocked country:
# But the files are read anyway:
all_EEZ = gpd.read_file('eez_v11.gpkg')
all_EEZ = gpd.read_file('eez_boundaries_v11.gpkg')
In [5]:
## !!! For presentation purposes !!! ###
# Reading the Slovakia Shape from data "country_shapes.geojson" file:
url = 'country_shapes (4).geojson'
countries = gpd.read_file(url).set_index("name")
#countries.tail(5)
gebco_svk = "GEBCO_2014_2D-SK.nc"
gebco = rasterio.open(gebco_svk)
band = gebco.read(1)
fig, ax = plt.subplots(figsize=(8, 12))
countries.loc[["SK"]].plot(ax=ax, color="none")
show(band, transform=gebco.transform, cmap="RdBu_r", ax=ax)
Out[5]:
<Axes: >
First, a land eligibility analysis, based on the criteria listede belo, is performed as follows:¶
Criteria for onshore wind implementation:
- 10km distance to airports
- 300m distance to major roads
- no natural protection areas
- maximum elevation of 2000m
- 1000m distance to built-up areas
- only on suitable land cover classes
Criteria for solar implementation:
only on suitable land cover classes
no natural protection areas
In [6]:
# Defining Exclusion and Inclusion areas as data frames:
Exclusion_area = pd.DataFrame({"land classes" :[111, 112, 121, 122, 123, 124, 142, 511, 512, 131, 322],
"code value" : [1, 2, 3, 4, 5, 6, 11, 40, 41, 7, 27],},
index = ['Continuous urban fabric', 'Discontinuous urban fabric',
'Industrial or commercial units', 'Road and rail networks and associated land',
'Port areas', 'Airports', 'Sport and leisure facilities', 'Water courses', 'Water bodies','Mineral extraction sites','Moors and heathland'])
Exclusion_area = Exclusion_area.rename_axis(index='descriptive names')
Inclusion_area = pd.DataFrame({"land classes" :[211, 212, 213, 231, 241, 242, 243, 321, 323, 324, 333],
"code value" : [12, 13, 14, 18, 19, 20, 21, 26, 28, 29, 32],},
index = ['Non-irrigated arable land', 'Permanently irrigated land', 'Rice fields',
'Pastures', 'Annual crops associated with permanent crops', 'Complex cultivation patterns',
'Land principally occupied by agriculture with significant areas of natural vegetation',
'Natural grasslands', 'Sclerophyllous vegetation', 'Transitional woodland-shrub',
'Sparsely vegetated areas'])
Inclusion_area = Inclusion_area.rename_axis(index='descriptive names')
# Defining function to plot the eligible area on the country map based on excluded & included areas:
def plot_area(masked, transform, shape):
fig, ax = plt.subplots(figsize=(5, 5))
ax = show(masked, transform=transform, cmap='Greens', vmin=0, ax=ax)
shape.plot(ax=ax, edgecolor='k', color='None', linewidth=2)
# Defining a function to calculate the eligible area portion of the country in percentage based on excluded and included areas:
def eligible_area(masked, excluder, shape):
return masked.sum() * float(excluder.res**2) / shape.geometry.area.sum() * 100
####################################################################################################
# Excluding the unsuitable area for onshore wind technology:
excluder_onshore_wind = ExclusionContainer()
excluder_onshore_wind.add_geometry("ne_10m_roads.gpkg", buffer=300) # 300m distance to major roads
excluder_onshore_wind.add_geometry("ne_10m_airports (1).gpkg", buffer=10000) # 10km distance to airports
excluder_onshore_wind.add_raster('WDPA_Oct2022_Public_shp-SVK.tif', crs=3035) # no natural protection areas
excluder_onshore_wind.add_raster('PROBAV_LC100_global_v3_0_1_2019_nrt_Discrete_Classification_map.tif', codes=[1,2,3,4,5,6,11,40,41], buffer=1000, crs=3035) #1000m distance to built-up areas
onshore_clases = [12, 13, 14, 18, 19, 20, 21, 26, 28, 29, 32] # only on suitable land cover classes
excluder_onshore_wind.add_raster('PROBAV_LC100_global_v3_0_1_2019_nrt_Discrete_Classification_map.tif', codes=onshore_clases, crs=3035)
masked, transform = shape_availability(shape_svk, excluder_onshore_wind)
# Calculating the eligible area for onshore wind technology in percentage (%):
onshore_eligible_area = eligible_area(masked, excluder_onshore_wind, shape_svk)
print()
print(f'Almost {round(onshore_eligible_area,3)} % of Slovakia is eligible for onshore wind technology development')
print()
# Plotting the eligible area for onshore wind technology:
plot_area(masked, transform, shape_svk)
plt.text(0.5, -0.4,f'Almost {round(onshore_eligible_area,3)} % of Slovakia is eligible for onshore wind technology development', ha='center', va='center', transform=plt.gca().transAxes)
plt.savefig('on-shore wind power eligibility.png')
plt.show()
Almost 5.432 % of Slovakia is eligible for onshore wind technology development
In [7]:
# Excluding the unsuitable area for solar power plants:
excluder_solar = ExclusionContainer()
Solar_classes = [1, 2, 3, 4, 5, 6, 11, 18, 26, 27, 30, 31, 41]
excluder_solar.add_raster('PROBAV_LC100_global_v3_0_1_2019_nrt_Discrete_Classification_map.tif', codes=Solar_classes, crs=3035)
excluder_solar = ExclusionContainer()
excluder_solar.add_raster('WDPA_Oct2022_Public_shp-SVK.tif', crs=3035)
masked, transform = shape_availability(shape_svk, excluder_solar)
# Calculating the eligible area for solar development in percentage (%):
eligible_area_solar = eligible_area(masked, excluder_solar,shape_svk)
print()
print(f'Almost {round(eligible_area_solar,2)} % of Slovakia land is eligible for solar')
print()
# Plotting the eligible area for solar development:
plot_area(masked, transform, shape_svk)
plt.text(0.5, -0.4,f'Almost {round(eligible_area_solar,2)} % of Slovakia land is eligible for solar implementation.', ha='center', va='center', transform=plt.gca().transAxes)
plt.savefig('Solar power eligibility.png')
plt.show()
Almost 62.69 % of Slovakia land is eligible for solar
Secondly, another land eligibility has been performed as follows, but this time using historical weather data of ERA5 dataset & atlite library.¶
The historical weather data of Slovakia in 2017 has been downloaded into an atlite.Cutout.
A buffer of 0.25 degreesdhas also been added to the geographical bounds ofthed countr:).
In [8]:
# Define the file path for the data, being downloaded from the ERA5 dataset:
data_file_path = "era5-2017-SVK.nc"
# defining the bound limits to create a spatial extent:
minx, miny, maxx, maxy = svk_regions.total_bounds
buffer = 0.25 # A buffer of 0.25 degrees has been considered as requested.
# Creating the cutout for the year 2017:
cutout = atlite.Cutout(
path=data_file_path,
module="era5",
x=slice(minx - buffer, maxx + buffer),
y=slice(miny - buffer, maxy + buffer),
time="2017",
)
# Calling the function cutout.prepare() & initiate the download:
cutout.prepare()
# the downloaded data is accessible here below :
cutout.data
INFO:atlite.data:Cutout already prepared.
Out[8]:
<xarray.Dataset>
Dimensions: (x: 25, y: 10, time: 8760)
Coordinates:
* x (x) float64 16.75 17.0 17.25 17.5 ... 22.25 22.5 22.75
* y (y) float64 47.5 47.75 48.0 48.25 ... 49.25 49.5 49.75
* time (time) datetime64[ns] 2017-01-01 ... 2017-12-31T23:00:00
lon (x) float64 dask.array<chunksize=(25,), meta=np.ndarray>
lat (y) float64 dask.array<chunksize=(10,), meta=np.ndarray>
Data variables: (12/13)
height (y, x) float32 dask.array<chunksize=(10, 25), meta=np.ndarray>
wnd100m (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
wnd_azimuth (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
roughness (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
influx_toa (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
influx_direct (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
... ...
albedo (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
solar_altitude (time, y, x) float64 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
solar_azimuth (time, y, x) float64 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
temperature (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
soil temperature (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
runoff (time, y, x) float32 dask.array<chunksize=(100, 10, 25), meta=np.ndarray>
Attributes:
module: era5
prepared_features: ['temperature', 'influx', 'wind', 'runoff', 'height']
chunksize_time: 100
Conventions: CF-1.6
history: 2024-01-31 15:41:18 GMT by grib_to_netcdf-2.25.1: /op...In [9]:
# Which data is contained in the cutout:
cutout
Out[9]:
<Cutout "era5-2017-SVK"> x = 16.75 ⟷ 22.75, dx = 0.25 y = 47.50 ⟷ 49.75, dy = 0.25 time = 2017-01-01 ⟷ 2017-12-31, dt = H module = era5 prepared_features = ['height', 'wind', 'influx', 'temperature', 'runoff']
In [10]:
# Included weather variables list:
cutout.prepared_features
Out[10]:
module feature
era5 height height
wind wnd100m
wind wnd_azimuth
wind roughness
influx influx_toa
influx influx_direct
influx influx_diffuse
influx albedo
influx solar_altitude
influx solar_azimuth
temperature temperature
temperature soil temperature
runoff runoff
dtype: object
In [11]:
# Creating availability matrixes for solar:
?cutout.availabilitymatrix
solar_AM = cutout.availabilitymatrix(shape_svk, excluder_solar)
#display(solar_AM)
# Creating availability matrixes for onshore wind:
onshore_AM = cutout.availabilitymatrix(shape_svk, excluder_onshore_wind)
#display(onshore_AM)
cap_per_sqkm = 3 # [MW/km2] For both wind and solar, a deployment density of 3 MW/km2 has been assumed.
# Calculating the area :
country_area = cutout.grid.set_index(["y", "x"]).to_crs(3035).area / 1e6
country_area = xr.DataArray(country_area, dims=("spatial"))
#display(country_area)
Compute availability matrix: 100%|███████████████████████████████████████████████| 8/8 [00:01<00:00, 4.56 gridcells/s] Compute availability matrix: 100%|███████████████████████████████████████████████| 8/8 [01:43<00:00, 12.97s/ gridcells]
Signature: cutout.availabilitymatrix( shapes, excluder, nprocesses=None, disable_progressbar=False, ) Docstring: Compute the eligible share within cutout cells in the overlap with shapes. For parallel calculation (nprocesses not None) the excluder must not be initialized and all raster references must be strings. Otherwise processes are colliding when reading from one common rasterio.DatasetReader. Parameters ---------- cutout : atlite.Cutout Cutout which the availability matrix is aligned to. shapes : geopandas.Series/geopandas.DataFrame Geometries for which the availabilities are calculated. excluder : atlite.gis.ExclusionContainer Container of all meta data or objects which to exclude, i.e. rasters and geometries. nprocesses : int, optional Number of processes to use for calculating the matrix. The paralle- lization can heavily boost the calculation speed. The default is None. disable_progressbar: bool, optional Disable the progressbar if nprocesses is not None. Then the `map` function instead of the `imap` function is used for the multiprocessing pool. This speeds up the calculation. Returns ------- availabilities : xr.DataArray DataArray of shape (|shapes|, |y|, |x|) containing all the eligible share of cutout cell (x,y) in the overlap with shape i. Notes ----- The rasterio (or GDAL) average downsampling returns different results dependent on how the target raster (the cutout raster) is spanned. Either it is spanned from the top left going downwards, e.g. Affine(0.25, 0, 0, 0, -0.25, 50), or starting in the lower left corner and going up, e.g. Affine(0.25, 0, 0, 0, 0.25, 50). Here we stick to the top down version which is why we use `cutout.transform_r` and flipping the y-axis in the end. File: c:\users\mahtab\anaconda3\envs\esm-2023\lib\site-packages\atlite\gis.py Type: method
In [12]:
# Creating capacity matrix for onshore wind:
capacity_matrix_onshore = onshore_AM.stack(spatial=["y", "x"]) * country_area * cap_per_sqkm
#display(capacity_matrix_onshore)
# Wind profile:
wind_onshore = cutout.wind(
atlite.windturbines.Vestas_V112_3MW,
matrix=capacity_matrix_onshore,
index= svk_regions.index,
per_unit=True,
)
svk_wind_onshore_2017_potential_time_series = wind_onshore.to_pandas()
display(svk_wind_onshore_2017_potential_time_series.head(15))
INFO:atlite.convert:Convert and aggregate 'wind'.
| NAME_1 | Banskobystrický | Bratislavský | Košický | Nitriansky | Prešovský | Trenčiansky | Trnavský | Žilinský |
|---|---|---|---|---|---|---|---|---|
| time | ||||||||
| 2017-01-01 00:00:00 | 0.000003 | 0.000060 | 0.007087 | 0.000000e+00 | 0.030756 | 0.000335 | 0.000001 | 0.013306 |
| 2017-01-01 01:00:00 | 0.000000 | 0.000005 | 0.007612 | 0.000000e+00 | 0.034854 | 0.000475 | 0.000040 | 0.014443 |
| 2017-01-01 02:00:00 | 0.000000 | 0.000066 | 0.007142 | 1.130307e-06 | 0.039974 | 0.000691 | 0.000331 | 0.015357 |
| 2017-01-01 03:00:00 | 0.000000 | 0.000171 | 0.003622 | 5.924064e-07 | 0.041654 | 0.001257 | 0.000336 | 0.017744 |
| 2017-01-01 04:00:00 | 0.000000 | 0.000144 | 0.003500 | 0.000000e+00 | 0.047368 | 0.001554 | 0.000341 | 0.020749 |
| 2017-01-01 05:00:00 | 0.000000 | 0.000156 | 0.003398 | 0.000000e+00 | 0.054182 | 0.002014 | 0.000318 | 0.027655 |
| 2017-01-01 06:00:00 | 0.000000 | 0.000167 | 0.003533 | 4.355234e-08 | 0.062805 | 0.002464 | 0.000342 | 0.035510 |
| 2017-01-01 07:00:00 | 0.000000 | 0.000406 | 0.004653 | 6.033786e-07 | 0.067074 | 0.002595 | 0.000470 | 0.042285 |
| 2017-01-01 08:00:00 | 0.000000 | 0.000687 | 0.004192 | 2.199573e-06 | 0.066282 | 0.002300 | 0.001231 | 0.044154 |
| 2017-01-01 09:00:00 | 0.000000 | 0.000343 | 0.005305 | 1.851060e-07 | 0.058913 | 0.001037 | 0.000373 | 0.044615 |
| 2017-01-01 10:00:00 | 0.000019 | 0.000218 | 0.005946 | 1.088100e-06 | 0.053568 | 0.000092 | 0.000304 | 0.027705 |
| 2017-01-01 11:00:00 | 0.000000 | 0.000181 | 0.007191 | 1.481370e-06 | 0.049395 | 0.000112 | 0.000274 | 0.025161 |
| 2017-01-01 12:00:00 | 0.000000 | 0.000299 | 0.008346 | 1.112223e-04 | 0.067957 | 0.000216 | 0.000336 | 0.031990 |
| 2017-01-01 13:00:00 | 0.000000 | 0.000523 | 0.008210 | 1.796761e-04 | 0.085420 | 0.000869 | 0.000486 | 0.044031 |
| 2017-01-01 14:00:00 | 0.000106 | 0.001866 | 0.008984 | 3.887791e-03 | 0.107708 | 0.005368 | 0.007772 | 0.060552 |
In [13]:
# Plotting the on-shore wind potential time series for Jan-2017:
plt.figure(figsize=(15, 6))
region_colors = {
'Banskobystrický': 'b',
'Bratislavský': 'g',
'Košický': 'r',
'Nitriansky': 'y',
'Prešovský': 'brown',
'Trenčiansky': 'pink',
'Trnavský': 'k',
'Žilinský': 'c',
}
for region, color in region_colors.items():
svk_wind_onshore_2017_potential_time_series.loc['2017-01', region].plot(color=color, label=region)
plt.ylabel('On-shore Wind Potential [p.u.]')
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] of all regions in the January, 2017')
plt.savefig('On-shore Wind Potential-January.png')
plt.legend()
plt.grid(True)
plt.show()
################################################
# Plotting the on-shore wind potential time series: ( Banskobystrický)
plt.figure(figsize=(12, 6))
svk_wind_onshore_2017_potential_time_series.iloc[:,0].plot(color='red')
plt.ylabel('On-shore Wind Potential [p.u.]')
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] for Banskobystrický region in the country of Slovokia')
plt.savefig('On-shore Wind Potential-Banskobystrický.png')
plt.show()
##################################################
# Plotting the on-shore wind potential time series: ( Bratislavský)
plt.figure(figsize=(12, 6))
svk_wind_onshore_2017_potential_time_series.iloc[:,1].plot(color='red')
plt.ylabel('On-shore Wind Potential [p.u.]')
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] for Bratislavský region in the country of Slovokia')
plt.savefig('On-shore Wind Potential-Bratislavský.png')
plt.show()
################################################
# Plotting the on-shore wind potential time series: ( Košický)
plt.figure(figsize=(12, 6))
svk_wind_onshore_2017_potential_time_series.iloc[:,2].plot(color='green')
plt.ylabel('On-shore Wind Potential [p.u.]')
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] for Košický region in the country of Slovokia')
plt.savefig('On-shore Wind Potential-Košický.png')
plt.show()
################################################
# Plotting the on-shore wind potential time series: ( Nitriansky)
plt.figure(figsize=(12, 6))
svk_wind_onshore_2017_potential_time_series.iloc[:,3].plot(color='grey')
plt.ylabel('On-shore Wind Potential [p.u.]')
plt.ylim(0, 1)
plt.title('On-shore Wind Potential [p.u.] for Nitriansky region in the country of Slovokia')
plt.savefig('On-shore Wind Potential-Nitriansky.png')
plt.show()
In [14]:
# Creating capacity matrix for solar:
capacity_matrix_solar = solar_AM.stack(spatial=['y', 'x']) * country_area * cap_per_sqkm
#display(capacity_matrix_solar)
# Solar PV profile:
pv = cutout.pv(
panel=atlite.solarpanels.CdTe,
matrix=capacity_matrix_solar,
orientation="latitude_optimal",
index= svk_regions.index,
per_unit=True,
)
svk_pv_2017_potential_time_series = pv.to_pandas()
svk_pv_2017_potential_time_series.tail(15)
INFO:atlite.convert:Convert and aggregate 'pv'.
[########################################] | 100% Completed | 5.36 s
Out[14]:
| NAME_1 | Banskobystrický | Bratislavský | Košický | Nitriansky | Prešovský | Trenčiansky | Trnavský | Žilinský |
|---|---|---|---|---|---|---|---|---|
| time | ||||||||
| 2017-12-31 09:00:00 | 0.028286 | 0.013923 | 0.074632 | 0.020029 | 0.142586 | 0.017215 | 0.014964 | 0.035910 |
| 2017-12-31 10:00:00 | 0.049915 | 0.022226 | 0.118199 | 0.042010 | 0.244426 | 0.039485 | 0.027253 | 0.065328 |
| 2017-12-31 11:00:00 | 0.156638 | 0.026005 | 0.158670 | 0.046129 | 0.304047 | 0.045949 | 0.020684 | 0.143634 |
| 2017-12-31 12:00:00 | 0.085264 | 0.050059 | 0.127709 | 0.042245 | 0.222590 | 0.027803 | 0.034511 | 0.067804 |
| 2017-12-31 13:00:00 | 0.050969 | 0.161593 | 0.068615 | 0.061506 | 0.163986 | 0.024984 | 0.110701 | 0.047229 |
| 2017-12-31 14:00:00 | 0.038435 | 0.237795 | 0.022509 | 0.117132 | 0.056720 | 0.030638 | 0.184676 | 0.016817 |
| 2017-12-31 15:00:00 | 0.012403 | 0.129751 | 0.005537 | 0.069491 | 0.010469 | 0.028245 | 0.108905 | 0.003444 |
| 2017-12-31 16:00:00 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 2017-12-31 17:00:00 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 2017-12-31 18:00:00 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 2017-12-31 19:00:00 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 2017-12-31 20:00:00 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 2017-12-31 21:00:00 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 2017-12-31 22:00:00 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 2017-12-31 23:00:00 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
In [15]:
# Plotting the solar potential time series for Jan-2017:
plt.figure(figsize=(15, 6))
for region, color in region_colors.items():
pv.loc['2017-01', region].plot(color=color, label=region)
plt.ylabel('On-shore Wind Potential [p.u.]')
plt.ylim(0, 1)
plt.title('Solar Potential [p.u.] of all regions in the January, 2017')
plt.savefig('Solar Potential-January.png')
plt.legend()
plt.grid(True)
plt.show()
##########################################################
# Plotting the solar potential time series: (only Banskobystrický)
plt.figure(figsize=(12, 6))
pv.to_pandas().iloc[:,0].plot(color='blue')
plt.ylabel('Solar Potential [p.u.]')
plt.ylim(0, 1)
plt.title('Solar Potential [p.u.] for Banskobystrický region in the country of Slovokia')
plt.savefig('Solar Potential-Banskobystrický.png')
plt.show()
##################################################
# Plotting the solar potential time series: (only Bratislavský)
plt.figure(figsize=(12, 6))
pv.to_pandas().iloc[:,1].plot(color='red')
plt.ylabel('Solar Potential [p.u.]')
plt.ylim(0, 1)
plt.title('Solar Potential [p.u.] for Bratislavský region in the country of Slovokia')
plt.savefig('Solar Potential-Bratislavský.png')
plt.show()
######################################################
# Plotting the solar potential time series: (only Košický)
plt.figure(figsize=(12, 6))
pv.to_pandas().iloc[:,2].plot(color='green')
plt.ylabel('Solar Potential [p.u.]')
plt.ylim(0, 1)
plt.title('Solar Potential [p.u.] for Košický region in the country of Slovokia')
plt.savefig('Solar Potential-Košický.png')
plt.show()
################################################################
# Plotting the solar potential time series: (only Nitriansky)
plt.figure(figsize=(12, 6))
pv.to_pandas().iloc[:,3].plot(color='g')
plt.ylabel('Solar Potential [p.u.]')
plt.ylim(0, 1)
plt.title('Solar Potential [p.u.] for Nitriansky region in the country of Slovokia')
plt.savefig('Solar Potential-Nitriansky.png')
plt.show()
In [48]:
from matplotlib.cm import ScalarMappable
## Building the Pypsa Model to minimize total annual system costs:
#display(gdf_svk.head(3))
n = pypsa.Network() # defining a network
start = '2017-01-01 00:00:00'
end = '2017-12-31 23:00:00'
timestamps = pd.date_range(start=start, end=end, freq='H')
n.set_snapshots(timestamps)
# Defining the Buses:
regions_points = svk_regions.representative_point()
for index, row in svk_regions.iterrows(): # A bus (Node) is defined for each region based on its geometry (latitude & Longitude)
bus_name = index
if bus_name in n.buses.index:
n.remove("Bus", bus_name)
region_point = regions_points.loc[index]
n.add("Bus", bus_name, x=region_point.x, y=region_point.y)
display(n.buses)
# !!!! FOR PRESENTATION: !!!
fig, ax = plt.subplots(figsize=(13, 9))
svk_regions.plot(ax=ax, facecolor='lightblue', edgecolor='black', linewidth=1.5)
regions_points.plot(ax=ax, color="black", markersize=20);
for bus_name, bus in n.buses.iterrows():
ax.annotate(bus_name, (bus['x'], bus['y']), textcoords="offset points", xytext=(3,10), ha='center')
plt.title('regions of Slovokia as Nodes in the Model')
plt.savefig('Nodes.png')
plt.show()
##################################################################
# !!!! FOR PRESENTATION: !!!
fig, ax = plt.subplots(figsize=(13, 9))
svk_regions.plot(ax=ax, facecolor='lightblue', edgecolor='black', linewidth=1.5)
regions_points.plot(ax=ax, color="lightblue", markersize=0.5);
for bus_name, bus in n.buses.iterrows():
ax.annotate(bus_name, (bus['x'], bus['y']), textcoords="offset points", xytext=(3,10), ha='center')
power_plants_scatter = gdf_svk.plot(
ax=ax ,
column = gdf_svk['primary_fuel'],
markersize = gdf_svk['capacity_mw'] / 1e1 *4,
legend= True,
cmap='viridis',
)
plt.title('Power plants Capacities per region in Slovokia')
plt.savefig('Power plants capacities.png')
plt.show()
| attribute | v_nom | type | x | y | carrier | unit | v_mag_pu_set | v_mag_pu_min | v_mag_pu_max | control | sub_network |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Bus | |||||||||||
| Banskobystrický | 1.0 | 19.424031 | 48.499321 | AC | 1.0 | 0.0 | inf | PQ | |||
| Bratislavský | 1.0 | 17.175805 | 48.331983 | AC | 1.0 | 0.0 | inf | PQ | |||
| Košický | 1.0 | 21.292701 | 48.676239 | AC | 1.0 | 0.0 | inf | PQ | |||
| Nitriansky | 1.0 | 18.336016 | 48.223509 | AC | 1.0 | 0.0 | inf | PQ | |||
| Prešovský | 1.0 | 21.219622 | 49.114946 | AC | 1.0 | 0.0 | inf | PQ | |||
| Trenčiansky | 1.0 | 18.215575 | 48.902964 | AC | 1.0 | 0.0 | inf | PQ | |||
| Trnavský | 1.0 | 17.645627 | 48.323494 | AC | 1.0 | 0.0 | inf | PQ | |||
| Žilinský | 1.0 | 19.196299 | 49.177544 | AC | 1.0 | 0.0 | inf | PQ |
In [17]:
year = 2040 # A year for cost projections has been selected
url = f"https://raw.githubusercontent.com/PyPSA/technology-data/master/outputs/costs_{year}.csv"
tech_properties = pd.read_csv(url)
#display(tech_properties.head(10))
tech_properties_svk = tech_properties[
(tech_properties['technology'] == 'nuclear') |
(tech_properties['technology'] == 'CCGT') |
(tech_properties['technology'] == 'coal') |
(tech_properties['technology'] == 'hydro') |
(tech_properties['technology'] == 'solar') |
(tech_properties['technology'] == 'gas') |
(tech_properties['technology'] == 'electrolysis') |
(tech_properties['technology'] == 'fuel cell') |
(tech_properties['technology'] == 'hydrogen storage underground') |
(tech_properties['technology'] == 'battery inverter') |
(tech_properties['technology'] == 'battery storage')
]
tech_properties_svk['technology'] = tech_properties_svk['technology'].replace('gas', 'CCGT')
display(tech_properties_svk)
| technology | parameter | value | unit | source | further description | |
|---|---|---|---|---|---|---|
| 29 | CCGT | FOM | 3.3006 | %/year | Danish Energy Agency, technology_data_for_el_a... | 05 Gas turb. CC, steam extract.: Fixed O&M |
| 30 | CCGT | VOM | 4.1000 | EUR/MWh | Danish Energy Agency, technology_data_for_el_a... | 05 Gas turb. CC, steam extract.: Variable O&M |
| 31 | CCGT | c_b | 2.1000 | 50oC/100oC | Danish Energy Agency, technology_data_for_el_a... | 05 Gas turb. CC, steam extract.: Cb coefficient |
| 32 | CCGT | c_v | 0.1500 | 50oC/100oC | Danish Energy Agency, technology_data_for_el_a... | 05 Gas turb. CC, steam extract.: Cv coefficient |
| 33 | CCGT | efficiency | 0.5900 | per unit | Danish Energy Agency, technology_data_for_el_a... | 05 Gas turb. CC, steam extract.: Electricity ... |
| 34 | CCGT | investment | 815.0000 | EUR/kW | Danish Energy Agency, technology_data_for_el_a... | 05 Gas turb. CC, steam extract.: Nominal inve... |
| 35 | CCGT | lifetime | 25.0000 | years | Danish Energy Agency, technology_data_for_el_a... | 05 Gas turb. CC, steam extract.: Technical li... |
| 412 | battery inverter | FOM | 0.5400 | %/year | Danish Energy Agency, technology_data_catalogu... | : Fixed O&M |
| 413 | battery inverter | efficiency | 0.9600 | per unit | Danish Energy Agency, technology_data_catalogu... | : Round trip efficiency DC |
| 414 | battery inverter | investment | 100.0000 | EUR/kW | Danish Energy Agency, technology_data_catalogu... | : Output capacity expansion cost investment |
| 415 | battery inverter | lifetime | 10.0000 | years | Danish Energy Agency, technology_data_catalogu... | : Technical lifetime |
| 416 | battery storage | investment | 94.0000 | EUR/kWh | Danish Energy Agency, technology_data_catalogu... | : Energy storage expansion cost investment |
| 417 | battery storage | lifetime | 30.0000 | years | Danish Energy Agency, technology_data_catalogu... | : Technical lifetime |
| 579 | coal | CO2 intensity | 0.3361 | tCO2/MWh_th | Entwicklung der spezifischen Kohlendioxid-Emis... | NaN |
| 580 | coal | FOM | 1.3100 | %/year | Lazard's levelized cost of energy analysis - v... | Calculated based on average of listed range, i... |
| 581 | coal | VOM | 3.3278 | EUR/MWh_e | Lazard's levelized cost of energy analysis - v... | Calculated based on average of listed range, i... |
| 582 | coal | efficiency | 0.3300 | p.u. | Lazard's levelized cost of energy analysis - v... | Calculated based on average of listed range, i... |
| 583 | coal | fuel | 9.2743 | EUR/MWh_th | DIW (2013): Current and propsective costs of e... | Based on IEA 2011 data, 99 USD/t. |
| 584 | coal | investment | 3905.3074 | EUR/kW_e | Lazard's levelized cost of energy analysis - v... | Higher costs include coal plants with CCS, but... |
| 585 | coal | lifetime | 40.0000 | years | Lazard's levelized cost of energy analysis - v... | NaN |
| 707 | electrolysis | FOM | 4.0000 | %/year | Danish Energy Agency, data_sheets_for_renewabl... | 86 AEC 100 MW: Fixed O&M |
| 708 | electrolysis | efficiency | 0.6532 | per unit | Danish Energy Agency, data_sheets_for_renewabl... | 86 AEC 100 MW: Hydrogen Output |
| 709 | electrolysis | efficiency-heat | 0.1849 | per unit | Danish Energy Agency, data_sheets_for_renewabl... | 86 AEC 100 MW: - hereof recoverable for dist... |
| 710 | electrolysis | investment | 384.9356 | EUR/kW_e | Danish Energy Agency, data_sheets_for_renewabl... | 86 AEC 100 MW: Specific investment |
| 711 | electrolysis | lifetime | 25.0000 | years | Danish Energy Agency, data_sheets_for_renewabl... | 86 AEC 100 MW: Technical lifetime |
| 712 | fuel cell | FOM | 5.0000 | %/year | Danish Energy Agency, technology_data_for_el_a... | 12 LT-PEMFC CHP: Fixed O&M |
| 713 | fuel cell | c_b | 1.2500 | 50oC/100oC | Danish Energy Agency, technology_data_for_el_a... | 12 LT-PEMFC CHP: Cb coefficient |
| 714 | fuel cell | efficiency | 0.5000 | per unit | Danish Energy Agency, technology_data_for_el_a... | 12 LT-PEMFC CHP: Electricity efficiency, annu... |
| 715 | fuel cell | investment | 950.0000 | EUR/kW_e | Danish Energy Agency, technology_data_for_el_a... | 12 LT-PEMFC CHP: Nominal investment |
| 716 | fuel cell | lifetime | 10.0000 | years | Danish Energy Agency, technology_data_for_el_a... | 12 LT-PEMFC CHP: Technical lifetime |
| 717 | CCGT | CO2 intensity | 0.1980 | tCO2/MWh_th | Stoichiometric calculation with 50 GJ/t CH4 | NaN |
| 718 | CCGT | fuel | 23.8481 | EUR/MWh_th | DIW (2013): Current and propsective costs of e... | Based on IEA 2011 data. |
| 745 | hydro | FOM | 1.0000 | %/year | DIW DataDoc http://hdl.handle.net/10419/80348 | from old pypsa cost assumptions |
| 746 | hydro | efficiency | 0.9000 | per unit | DIW DataDoc http://hdl.handle.net/10419/80348 | from old pypsa cost assumptions |
| 747 | hydro | investment | 2208.1616 | EUR/kWel | DIW DataDoc http://hdl.handle.net/10419/80348 | from old pypsa cost assumptions |
| 748 | hydro | lifetime | 80.0000 | years | IEA2010 | from old pypsa cost assumptions |
| 760 | hydrogen storage underground | FOM | 0.0000 | %/year | Danish Energy Agency, technology_data_catalogu... | 151c Hydrogen Storage - Caverns: Fixed O&M |
| 761 | hydrogen storage underground | VOM | 0.0000 | EUR/MWh | Danish Energy Agency, technology_data_catalogu... | 151c Hydrogen Storage - Caverns: Variable O&M |
| 762 | hydrogen storage underground | investment | 1.5000 | EUR/kWh | Danish Energy Agency, technology_data_catalogu... | 151c Hydrogen Storage - Caverns: Specific inv... |
| 763 | hydrogen storage underground | lifetime | 100.0000 | years | Danish Energy Agency, technology_data_catalogu... | 151c Hydrogen Storage - Caverns: Technical li... |
| 814 | nuclear | FOM | 1.2700 | %/year | Lazard's levelized cost of energy analysis - v... | U.S. specific costs including newly commission... |
| 815 | nuclear | VOM | 3.6188 | EUR/MWh_e | Lazard's levelized cost of energy analysis - v... | U.S. specific costs including newly commission... |
| 816 | nuclear | efficiency | 0.3260 | p.u. | Lazard's levelized cost of energy analysis - v... | Based on heat rate of 10.45 MMBtu/MWh_e and 3.... |
| 817 | nuclear | fuel | 3.3122 | EUR/MWh_th | DIW (2013): Current and propsective costs of e... | Based on IEA 2011 data. |
| 818 | nuclear | investment | 8769.6136 | EUR/kW_e | Lazard's levelized cost of energy analysis - v... | U.S. specific costs including newly commission... |
| 819 | nuclear | lifetime | 40.0000 | years | Lazard's levelized cost of energy analysis - v... | NaN |
| 858 | solar | FOM | 2.0400 | %/year | Calculated. See 'further description'. | Mixed investment costs based on average of 50%... |
| 859 | solar | VOM | 0.0100 | EUR/MWhel | RES costs made up to fix curtailment order | from old pypsa cost assumptions |
| 860 | solar | investment | 407.8706 | EUR/kW_e | Calculated. See 'further description'. | Mixed investment costs based on average of 50%... |
| 861 | solar | lifetime | 40.0000 | years | Calculated. See 'further description'. | Mixed investment costs based on average of 50%... |
In [18]:
# Dataframe of the technologies:
df_CCGT = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'CCGT']
df_CCGT.set_index('parameter', inplace=True)
#display(df_CCGT)
df_solar = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'solar']
df_solar.set_index('parameter', inplace=True)
df_solar.loc['fuel'] = ['solar', 0, None, None, None]
df_solar.loc['efficiency'] = ['solar', 1, None, None, None]
#display(df_solar)
df_nuclear = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'nuclear']
df_nuclear.set_index('parameter', inplace=True)
#display(df_nuclear)
df_hydro = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'hydro']
df_hydro.set_index('parameter', inplace=True)
df_hydro.loc['fuel'] = ['hydro', 0, None, None, None]
df_hydro.loc['VOM'] = ['hydro', 0, None, None, None]
#display(df_hydro)
df_coal = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'coal']
df_coal.set_index('parameter', inplace=True)
#display(df_coal)
df_electrolysis = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'electrolysis']
df_electrolysis.set_index('parameter', inplace=True)
df_fuel_cell = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'fuel cell']
df_fuel_cell.set_index('parameter', inplace=True)
df_hydrogen_storage_underground = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'hydrogen storage underground']
df_hydrogen_storage_underground.set_index('parameter', inplace=True)
display(df_electrolysis)
df_battery_inverter = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'battery inverter']
df_battery_inverter.set_index('parameter', inplace=True)
df_battery_storage = tech_properties_svk[tech_properties_svk.iloc[:, 0] == 'battery storage']
df_battery_storage.set_index('parameter', inplace=True)
| technology | value | unit | source | further description | |
|---|---|---|---|---|---|
| parameter | |||||
| FOM | electrolysis | 4.0000 | %/year | Danish Energy Agency, data_sheets_for_renewabl... | 86 AEC 100 MW: Fixed O&M |
| efficiency | electrolysis | 0.6532 | per unit | Danish Energy Agency, data_sheets_for_renewabl... | 86 AEC 100 MW: Hydrogen Output |
| efficiency-heat | electrolysis | 0.1849 | per unit | Danish Energy Agency, data_sheets_for_renewabl... | 86 AEC 100 MW: - hereof recoverable for dist... |
| investment | electrolysis | 384.9356 | EUR/kW_e | Danish Energy Agency, data_sheets_for_renewabl... | 86 AEC 100 MW: Specific investment |
| lifetime | electrolysis | 25.0000 | years | Danish Energy Agency, data_sheets_for_renewabl... | 86 AEC 100 MW: Technical lifetime |
In [19]:
# Calculating the marginal costs of the technologies:
marginalcost_Coal = df_coal.loc['fuel', 'value'] / df_coal.loc['efficiency', 'value'] + df_coal.loc['VOM', 'value']
marginalcost_Solar = df_solar.loc['fuel', 'value'] / df_solar.loc['efficiency', 'value'] + df_solar.loc['VOM', 'value']
marginalcost_Nuclear = df_nuclear.loc['fuel', 'value'] / df_nuclear.loc['efficiency', 'value'] + df_nuclear.loc['VOM', 'value']
marginalcost_Gas = df_CCGT.loc['fuel', 'value'] / df_CCGT.loc['efficiency', 'value'] + df_CCGT.loc['VOM', 'value']
marginalcost_Hydro = df_hydro.loc['fuel', 'value'] / df_hydro.loc['efficiency', 'value'] + df_hydro.loc['VOM', 'value']
marginalcost_Wind = 0
marginal_costs ={ 'Coal' : marginalcost_Coal,
'Solar' : marginalcost_Solar,
'Nuclear': marginalcost_Nuclear,
'Gas': marginalcost_Gas,
'Hydro' : marginalcost_Hydro,
'Wind' : marginalcost_Wind,
}
In [20]:
# Calculating the capital costs of the technologies:
r = 0.07
annuity_coal = (r / (1 - (1 + r) ** -df_coal.loc['lifetime', 'value'])) * (df_coal.loc['investment', 'value']) * 1000
annuity_solar = (r / (1 - (1 + r) ** -df_solar.loc['lifetime', 'value'])) * (df_solar.loc['investment', 'value']) * 1000
annuity_hydro = (r / (1 - (1 + r) ** -df_hydro.loc['lifetime', 'value'])) * (df_hydro.loc['investment', 'value']) * 1000
annuity_nuclear = (r / (1 - (1 + r) ** -df_nuclear.loc['lifetime', 'value'])) * (df_nuclear.loc['investment', 'value']) * 1000
annuity_CCGT = (r / (1 - (1 + r) ** -df_CCGT.loc['lifetime', 'value'])) * (df_CCGT.loc['investment', 'value']) * 1000
annuity_electrolysis = (r / (1 - (1 + r) ** -df_electrolysis.loc['lifetime', 'value'])) * (df_electrolysis.loc['investment', 'value']) * 1000
annuity_fuel_cell = (r / (1 - (1 + r) ** -df_fuel_cell.loc['lifetime', 'value'])) * (df_fuel_cell.loc['investment', 'value']) * 1000
annuity_hydrogen_storage_underground = (r / (1 - (1 + r) ** -df_hydrogen_storage_underground.loc['lifetime', 'value'])) * (df_hydrogen_storage_underground.loc['investment', 'value']) * 1000
annuity_battery_inverter = (r / (1 - (1 + r) ** -df_battery_inverter.loc['lifetime', 'value'])) * (df_battery_inverter.loc['investment', 'value']) * 1000
annuity_battery_storage = (r / (1 - (1 + r) ** -df_battery_storage.loc['lifetime', 'value'])) * (df_battery_storage.loc['investment', 'value']) * 1000
FOM_coal = (df_coal.loc['FOM', 'value'] / 100) * df_coal.loc['investment', 'value'] *1000
FOM_solar = (df_solar.loc['FOM', 'value'] / 100) * df_solar.loc['investment', 'value'] *1000
FOM_nuclear = (df_nuclear.loc['FOM', 'value'] / 100) * df_nuclear.loc['investment', 'value'] *1000
FOM_hydro = (df_hydro.loc['FOM', 'value'] / 100) * df_hydro.loc['investment', 'value'] *1000
FOM_CCGT = (df_CCGT.loc['FOM', 'value'] / 100) * df_CCGT.loc['investment', 'value'] *1000
FOM_electrolysis = (df_electrolysis.loc['FOM', 'value'] / 100) * df_electrolysis.loc['investment', 'value'] *1000
FOM_fuel_cell = (df_fuel_cell.loc['FOM', 'value'] / 100) * df_fuel_cell.loc['investment', 'value'] *1000
FOM_battery_inverter = (df_battery_inverter.loc['FOM', 'value'] / 100) * df_battery_inverter.loc['investment', 'value'] *1000
capitalcost_coal = annuity_coal + FOM_coal
capitalcost_solar = annuity_solar + FOM_solar
capitalcost_nuclear = annuity_nuclear + FOM_nuclear
capitalcost_hydro = annuity_hydro + FOM_hydro
capitalcost_CCGT = annuity_CCGT + FOM_CCGT
capitalcost_electrolysis= annuity_electrolysis + FOM_electrolysis
capitalcost_fuel_cell= annuity_fuel_cell + FOM_fuel_cell
capitalcost_hydrogen_storage_underground = annuity_hydrogen_storage_underground
capitalcost_battery_storage = annuity_battery_storage
capitalcost_battery_inverter= annuity_battery_inverter + FOM_battery_inverter
display(capitalcost_battery_storage)
7575.121930044452
In [21]:
## Defining the Generators:
power_plants_df_with_regions = gpd.sjoin(gdf_svk, svk_regions, how="left", op="within") # spatial join to specify in which region each power plant is located
#display(power_plants_df_with_regions.head(5))
aggregated_power_plants = power_plants_df_with_regions.groupby(['primary_fuel', 'index_right']).first() # grouping the power plant data based on technology and region
#display(aggregated_power_plants)
# Adding the Carriers & Corresponding emissions:
emissions = {
'Coal': 0.34,
'Gas': 0.2,
'Hydro': 0,
'Nuclear': 0,
'Solar': 0,
'Wind': 0
}
n.madd(
"Carrier",
["Coal", "Gas", "Hydro", "Wind", "Solar", "Nuclear"],
co2_emissions=emissions,
)
# Adding the generators to the network
for index, row in aggregated_power_plants.iterrows():
generator_name = f"{index[0]}_{index[1]}" # Specifying a name for each generator
if generator_name in n.generators.index: # Removing already created generators to avoid errors regarding duplication of generators name
n.remove("Generator", generator_name)
p_nom = row['capacity_mw']
if index[0]=='Hydro':
p_max_pu = (row['estimated_generation_gwh_2017'] * 1000) / (row['capacity_mw'] * 365 * 24 )
else:
p_max_pu = 1
n.add("Generator", # Adding generators
generator_name,
bus=index[1],
p_nom_extendable = False, # The Generators are not extendable
capital_cost = 0, # The Generators have no capital costs
p_nom=p_nom,
marginal_cost = marginal_costs[index[0]],
p_max_pu = p_max_pu,
carrier=index[0],
)
for generator_name in n.generators.index: # Disregarding existing wind and solar capacities
if "Wind" in generator_name or "Solar" in generator_name:
n.remove("Generator", generator_name)
display(n.generators.carrier)
Generator Coal_Košický Coal Coal_Trenčiansky Coal Gas_Bratislavský Gas Gas_Nitriansky Gas Gas_Trnavský Gas Hydro_Banskobystrický Hydro Hydro_Bratislavský Hydro Hydro_Košický Hydro Hydro_Trenčiansky Hydro Hydro_Trnavský Hydro Hydro_Žilinský Hydro Nuclear_Nitriansky Nuclear Nuclear_Trnavský Nuclear Name: carrier, dtype: object
In [22]:
load_time_series = pd.read_csv('load.csv', index_col=0, parse_dates=True) # Reading the load file
load_data_series_svk = load_time_series[['SK']].copy() # Cutting up the Slovakia part of the time series
#display(load_data_series_svk.head(10))
new_index = load_data_series_svk.index.strftime('%Y-%m-%d %H:%M:%S').str.replace('2013', '2017')
load_data_series_svk.index = pd.to_datetime(new_index)
load_data_series_svk = load_data_series_svk.rename_axis('Time')
#display(load_data_series_svk)
# Defining a dictionary including regions population [ref.: https://www.citypopulation.de/en/slovakia/cities/]:
cities_population = {'Banskobystrický':617777,
'Bratislavský' : 728370,
'Košický': 779505,
'Nitriansky': 670696,
'Prešovský': 808090,
'Trenčiansky': 570675,
'Trnavský': 565573,
'Žilinský': 688106}
total_population_skv = sum(cities_population.values()) # Calculating the total population of Slovokia
# Adding extra columns as region's load value time series to the primary load time series:
for region in cities_population.keys():
load_data_series_svk ['skv_load_per_capita'] = load_data_series_svk['SK'] / total_population_skv # calculating the load per capita
load_data_series_svk [f'load_{region}'] = load_data_series_svk ['skv_load_per_capita'] * cities_population [region] # Calculating the load on each region based on the region's population
#display(load_data_series_svk)
for index , bus in n.buses.iterrows():
Load_name = f'Demand_load_{index}'
if (n.loads.index == Load_name).any():
n.remove("Load", Load_name)
p_set_data = load_data_series_svk[f'load_{index}']
n.add("Load", # Adding Electricity loads for each bus (region)
Load_name,
bus=index,
p_set=p_set_data,
carrier="electricity",
)
display(n.loads_t.p_set)
ax = n.loads_t.p_set.iloc[:24].plot(figsize=(10, 6),ylim=[0, 800], ylabel="Load [MW]", xlabel="Time", title='Electricity Demand Load during a single day in different regions of Slovokia (01.01.2017)')
plt.savefig('electricity_demand_load.png')
plt.show()
#n.snapshots
| Load | Demand_load_Banskobystrický | Demand_load_Bratislavský | Demand_load_Košický | Demand_load_Nitriansky | Demand_load_Prešovský | Demand_load_Trenčiansky | Demand_load_Trnavský | Demand_load_Žilinský |
|---|---|---|---|---|---|---|---|---|
| snapshot | ||||||||
| 2017-01-01 00:00:00 | 364.359299 | 429.586052 | 459.745013 | 395.570447 | 476.604188 | 336.578964 | 333.569851 | 405.838708 |
| 2017-01-01 01:00:00 | 363.843475 | 428.977887 | 459.094152 | 395.010438 | 475.929459 | 336.102469 | 333.097616 | 405.264162 |
| 2017-01-01 02:00:00 | 362.529793 | 427.429033 | 457.436561 | 393.584226 | 474.211084 | 334.888948 | 331.894944 | 403.800928 |
| 2017-01-01 03:00:00 | 365.231299 | 430.614155 | 460.845294 | 396.517143 | 477.744817 | 337.384479 | 334.368164 | 406.809978 |
| 2017-01-01 04:00:00 | 381.691731 | 450.021296 | 481.614908 | 414.387582 | 499.276067 | 352.589897 | 349.437641 | 425.144300 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2017-12-31 19:00:00 | 460.693214 | 543.165440 | 581.298209 | 500.156360 | 602.614826 | 425.567964 | 421.763263 | 513.139474 |
| 2017-12-31 20:00:00 | 435.896276 | 513.929412 | 550.009674 | 473.235307 | 570.178918 | 402.661652 | 399.061740 | 485.519601 |
| 2017-12-31 21:00:00 | 405.975189 | 478.651922 | 512.255539 | 440.751170 | 531.040312 | 375.021879 | 371.669074 | 452.192237 |
| 2017-12-31 22:00:00 | 377.655121 | 445.262061 | 476.521552 | 410.005194 | 493.995935 | 348.861055 | 345.742136 | 420.648154 |
| 2017-12-31 23:00:00 | 372.171151 | 438.796364 | 469.601932 | 404.051465 | 486.822567 | 343.795207 | 340.721578 | 414.539877 |
8760 rows × 8 columns
In [23]:
# Defining regions area as a dictionary[km^2]:
Regions_area = {
'Banskobystrický' : 9455,
'Bratislavský' : 2053,
'Košický': 6753,
'Nitriansky': 6343,
'Prešovský' : 8974 ,
'Trenčiansky': 4502 ,
'Trnavský': 4145,
'Žilinský': 6809,
}
# Adding solar and on-shore wind generators per region:
n.buses
for index , bus in n.buses.iterrows():
solar_generator_name = f'Solar_{index}'
if solar_generator_name in n.generators.index:
n.remove("Generator", solar_generator_name)
n.add("Generator", # Adding solar generators
solar_generator_name,
bus=index,
p_nom_extendable = False, # The Generators are not extendable
capital_cost = 0, # The Generators have no capital costs
marginal_cost = marginalcost_Solar,
p_nom= Regions_area[index] * cap_per_sqkm,
p_max_pu = svk_pv_2017_potential_time_series[index],
carrier='Solar')
wind_generator_name = f'Wind_{index}'
if wind_generator_name in n.generators.index:
n.remove("Generator", wind_generator_name)
n.add("Generator", # Adding wind generators
wind_generator_name,
bus=index,
p_nom_extendable = False, # The Generators are not extendable
capital_cost = 0, # The Generators have no capital costs
marginal_cost = marginalcost_Wind,
p_nom= Regions_area[index] * cap_per_sqkm,
p_max_pu = svk_wind_onshore_2017_potential_time_series[index],
carrier='Wind',
)
display( n.generators.carrier.map(n.carriers.co2_emissions))
display(n.generators_t.p_max_pu)
Generator Coal_Košický 0.34 Coal_Trenčiansky 0.34 Gas_Bratislavský 0.20 Gas_Nitriansky 0.20 Gas_Trnavský 0.20 Hydro_Banskobystrický 0.00 Hydro_Bratislavský 0.00 Hydro_Košický 0.00 Hydro_Trenčiansky 0.00 Hydro_Trnavský 0.00 Hydro_Žilinský 0.00 Nuclear_Nitriansky 0.00 Nuclear_Trnavský 0.00 Solar_Banskobystrický 0.00 Wind_Banskobystrický 0.00 Solar_Bratislavský 0.00 Wind_Bratislavský 0.00 Solar_Košický 0.00 Wind_Košický 0.00 Solar_Nitriansky 0.00 Wind_Nitriansky 0.00 Solar_Prešovský 0.00 Wind_Prešovský 0.00 Solar_Trenčiansky 0.00 Wind_Trenčiansky 0.00 Solar_Trnavský 0.00 Wind_Trnavský 0.00 Solar_Žilinský 0.00 Wind_Žilinský 0.00 Name: carrier, dtype: float64
| Generator | Solar_Banskobystrický | Wind_Banskobystrický | Solar_Bratislavský | Wind_Bratislavský | Solar_Košický | Wind_Košický | Solar_Nitriansky | Wind_Nitriansky | Solar_Prešovský | Wind_Prešovský | Solar_Trenčiansky | Wind_Trenčiansky | Solar_Trnavský | Wind_Trnavský | Solar_Žilinský | Wind_Žilinský |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| snapshot | ||||||||||||||||
| 2017-01-01 00:00:00 | 0.0 | 0.000003 | 0.0 | 0.000060 | 0.0 | 0.007087 | 0.0 | 0.000000e+00 | 0.0 | 0.030756 | 0.0 | 0.000335 | 0.0 | 0.000001 | 0.0 | 0.013306 |
| 2017-01-01 01:00:00 | 0.0 | 0.000000 | 0.0 | 0.000005 | 0.0 | 0.007612 | 0.0 | 0.000000e+00 | 0.0 | 0.034854 | 0.0 | 0.000475 | 0.0 | 0.000040 | 0.0 | 0.014443 |
| 2017-01-01 02:00:00 | 0.0 | 0.000000 | 0.0 | 0.000066 | 0.0 | 0.007142 | 0.0 | 1.130307e-06 | 0.0 | 0.039974 | 0.0 | 0.000691 | 0.0 | 0.000331 | 0.0 | 0.015357 |
| 2017-01-01 03:00:00 | 0.0 | 0.000000 | 0.0 | 0.000171 | 0.0 | 0.003622 | 0.0 | 5.924064e-07 | 0.0 | 0.041654 | 0.0 | 0.001257 | 0.0 | 0.000336 | 0.0 | 0.017744 |
| 2017-01-01 04:00:00 | 0.0 | 0.000000 | 0.0 | 0.000144 | 0.0 | 0.003500 | 0.0 | 0.000000e+00 | 0.0 | 0.047368 | 0.0 | 0.001554 | 0.0 | 0.000341 | 0.0 | 0.020749 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2017-12-31 19:00:00 | 0.0 | 0.009669 | 0.0 | 0.135380 | 0.0 | 0.021920 | 0.0 | 9.197736e-02 | 0.0 | 0.127549 | 0.0 | 0.094018 | 0.0 | 0.114936 | 0.0 | 0.116168 |
| 2017-12-31 20:00:00 | 0.0 | 0.006298 | 0.0 | 0.135983 | 0.0 | 0.022566 | 0.0 | 5.302415e-02 | 0.0 | 0.130486 | 0.0 | 0.072025 | 0.0 | 0.096996 | 0.0 | 0.096775 |
| 2017-12-31 21:00:00 | 0.0 | 0.002918 | 0.0 | 0.151189 | 0.0 | 0.019083 | 0.0 | 5.068621e-02 | 0.0 | 0.128027 | 0.0 | 0.053257 | 0.0 | 0.105005 | 0.0 | 0.081824 |
| 2017-12-31 22:00:00 | 0.0 | 0.000884 | 0.0 | 0.114648 | 0.0 | 0.016145 | 0.0 | 4.555578e-02 | 0.0 | 0.118321 | 0.0 | 0.033475 | 0.0 | 0.077719 | 0.0 | 0.082753 |
| 2017-12-31 23:00:00 | 0.0 | 0.000821 | 0.0 | 0.158193 | 0.0 | 0.023487 | 0.0 | 4.725993e-02 | 0.0 | 0.119073 | 0.0 | 0.027659 | 0.0 | 0.121956 | 0.0 | 0.073361 |
8760 rows × 16 columns
In [24]:
# Comparing solar and wind potential in January 2017:
n.generators_t.p_max_pu.groupby(n.generators.carrier, axis=1).mean().iloc[:744].plot(ylabel="p.u.", figsize=(15, 6), color=['y', 'b'])
plt.savefig('Mean renewable energy capacity of implemented Generators-Jan-2017.png')
plt.show()
# Comparing solar and wind potential in 1.January 2017 (as a single day):
n.generators_t.p_max_pu.groupby(n.generators.carrier, axis=1).mean().iloc[:24].plot(ylabel="p.u.", figsize=(15, 6), color=['y', 'b'])
plt.savefig('Mean renewable energy capacity of implemented Generators-one day-1.Jan-2017.png')
plt.show()
In [25]:
# Difining a list of tuples, including neighboring regions to create transmission lines between them:
neighbor_buses = [('Bratislavský','Trnavský'),
('Trnavský','Nitriansky'),
('Trnavský','Trenčiansky'),
('Nitriansky','Trenčiansky'),
('Nitriansky','Banskobystrický'),
('Banskobystrický','Trenčiansky'),
('Žilinský','Trenčiansky'),
('Banskobystrický','Žilinský'),
('Prešovský','Žilinský'),
('Prešovský','Banskobystrický'),
('Košický','Banskobystrický'),
('Košický','Prešovský')]
for region1, region2 in neighbor_buses:
distance = np.sqrt((regions_points[region1].x - regions_points[region2].x) ** 2 + # Calculating the distance for each transmission line based on the representative points of buses
(regions_points[region1].y - regions_points[region2].y) ** 2)
length = 1.5 * distance # As requested, the length of 1.5 times the crow-fly distance between the regions’ representative points
cost = 400 * length # Assumption for costs of 400€/MW/km
# Adding bidirectional links between neighboring regions
TL1_name = f"TL1_{region1}_{region2}"
if TL1_name in n.links.index:
n.remove("Link", TL1_name)
n.add("Link",
TL1_name,
bus0=region1,
bus1=region2,
capital_cost=cost,
p_max_pu=1.0, # no transmission losses
efficiency=1.0, # no transmission losses
length=length,
)
TL2_name = f"TL2_{region2}_{region1}"
if TL2_name in n.links.index:
n.remove("Link", TL2_name)
n.add("Link",
TL2_name,
bus0=region2,
bus1=region1,
capital_cost=cost,
p_max_pu=1.0, # no transmission losses
efficiency=1.0, # no transmission losses
length=length,
)
display(n.links)
| attribute | bus0 | bus1 | type | carrier | efficiency | build_year | lifetime | p_nom | p_nom_extendable | p_nom_min | ... | shut_down_cost | min_up_time | min_down_time | up_time_before | down_time_before | ramp_limit_up | ramp_limit_down | ramp_limit_start_up | ramp_limit_shut_down | p_nom_opt |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Link | |||||||||||||||||||||
| TL1_Bratislavský_Trnavský | Bratislavský | Trnavský | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL2_Trnavský_Bratislavský | Trnavský | Bratislavský | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL1_Trnavský_Nitriansky | Trnavský | Nitriansky | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL2_Nitriansky_Trnavský | Nitriansky | Trnavský | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL1_Trnavský_Trenčiansky | Trnavský | Trenčiansky | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL2_Trenčiansky_Trnavský | Trenčiansky | Trnavský | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL1_Nitriansky_Trenčiansky | Nitriansky | Trenčiansky | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL2_Trenčiansky_Nitriansky | Trenčiansky | Nitriansky | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL1_Nitriansky_Banskobystrický | Nitriansky | Banskobystrický | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL2_Banskobystrický_Nitriansky | Banskobystrický | Nitriansky | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL1_Banskobystrický_Trenčiansky | Banskobystrický | Trenčiansky | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL2_Trenčiansky_Banskobystrický | Trenčiansky | Banskobystrický | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL1_Žilinský_Trenčiansky | Žilinský | Trenčiansky | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL2_Trenčiansky_Žilinský | Trenčiansky | Žilinský | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL1_Banskobystrický_Žilinský | Banskobystrický | Žilinský | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL2_Žilinský_Banskobystrický | Žilinský | Banskobystrický | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL1_Prešovský_Žilinský | Prešovský | Žilinský | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL2_Žilinský_Prešovský | Žilinský | Prešovský | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL1_Prešovský_Banskobystrický | Prešovský | Banskobystrický | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL2_Banskobystrický_Prešovský | Banskobystrický | Prešovský | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL1_Košický_Banskobystrický | Košický | Banskobystrický | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL2_Banskobystrický_Košický | Banskobystrický | Košický | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL1_Košický_Prešovský | Košický | Prešovský | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 | ||
| TL2_Prešovský_Košický | Prešovský | Košický | 1.0 | 0 | inf | 0.0 | False | 0.0 | ... | 0.0 | 0 | 0 | 1 | 0 | NaN | NaN | 1.0 | 1.0 | 0.0 |
24 rows × 32 columns
In [26]:
# Adding Battery storage units:
for bus in n.buses.index:
battery_storage_U_name = f"Battery_storage_{bus}"
if battery_storage_U_name in n.storage_units.index:
n.remove("StorageUnit", battery_storage_U_name)
n.add(
"StorageUnit",
battery_storage_U_name,
max_hours= 6,
capital_cost = 6 * capitalcost_battery_storage + capitalcost_battery_inverter,
bus = bus,
marginal_cost =0,
carrier = "battery storage",
efficiency_dispatch = df_battery_inverter.loc['efficiency', 'value'],
efficiency_store = df_battery_inverter.loc['efficiency', 'value'],
cyclic_state_of_charge = True,
p_nom_extendable = True,
)
n.storage_units
Out[26]:
| attribute | bus | control | type | p_nom | p_nom_extendable | p_nom_min | p_nom_max | p_min_pu | p_max_pu | p_set | ... | state_of_charge_initial_per_period | state_of_charge_set | cyclic_state_of_charge | cyclic_state_of_charge_per_period | max_hours | efficiency_store | efficiency_dispatch | standing_loss | inflow | p_nom_opt |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| StorageUnit | |||||||||||||||||||||
| Battery_storage_Banskobystrický | Banskobystrický | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.96 | 0.96 | 0.0 | 0.0 | 0.0 | |
| Battery_storage_Bratislavský | Bratislavský | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.96 | 0.96 | 0.0 | 0.0 | 0.0 | |
| Battery_storage_Košický | Košický | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.96 | 0.96 | 0.0 | 0.0 | 0.0 | |
| Battery_storage_Nitriansky | Nitriansky | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.96 | 0.96 | 0.0 | 0.0 | 0.0 | |
| Battery_storage_Prešovský | Prešovský | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.96 | 0.96 | 0.0 | 0.0 | 0.0 | |
| Battery_storage_Trenčiansky | Trenčiansky | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.96 | 0.96 | 0.0 | 0.0 | 0.0 | |
| Battery_storage_Trnavský | Trnavský | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.96 | 0.96 | 0.0 | 0.0 | 0.0 | |
| Battery_storage_Žilinský | Žilinský | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.96 | 0.96 | 0.0 | 0.0 | 0.0 |
8 rows × 29 columns
In [27]:
# Adding hydrogen storage units
for bus in n.buses.index:
storage_U_name = f"hydrogen_storage_{bus}"
if storage_U_name in n.storage_units.index:
n.remove("StorageUnit", storage_U_name)
n.add(
"StorageUnit",
storage_U_name,
max_hours= 336,
capital_cost = 336 * capitalcost_hydrogen_storage_underground+capitalcost_fuel_cell+capitalcost_electrolysis,
bus = bus,
marginal_cost = 0,
carrier = "hydrogen storage underground",
efficiency_dispatch = df_fuel_cell.loc['efficiency', 'value'],
efficiency_store = df_electrolysis.loc['efficiency', 'value'],
cyclic_state_of_charge = True,
p_nom_extendable = True,
)
display(n.storage_units)
| attribute | bus | control | type | p_nom | p_nom_extendable | p_nom_min | p_nom_max | p_min_pu | p_max_pu | p_set | ... | state_of_charge_initial_per_period | state_of_charge_set | cyclic_state_of_charge | cyclic_state_of_charge_per_period | max_hours | efficiency_store | efficiency_dispatch | standing_loss | inflow | p_nom_opt |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| StorageUnit | |||||||||||||||||||||
| Battery_storage_Banskobystrický | Banskobystrický | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.9600 | 0.96 | 0.0 | 0.0 | 0.0 | |
| Battery_storage_Bratislavský | Bratislavský | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.9600 | 0.96 | 0.0 | 0.0 | 0.0 | |
| Battery_storage_Košický | Košický | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.9600 | 0.96 | 0.0 | 0.0 | 0.0 | |
| Battery_storage_Nitriansky | Nitriansky | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.9600 | 0.96 | 0.0 | 0.0 | 0.0 | |
| Battery_storage_Prešovský | Prešovský | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.9600 | 0.96 | 0.0 | 0.0 | 0.0 | |
| Battery_storage_Trenčiansky | Trenčiansky | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.9600 | 0.96 | 0.0 | 0.0 | 0.0 | |
| Battery_storage_Trnavský | Trnavský | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.9600 | 0.96 | 0.0 | 0.0 | 0.0 | |
| Battery_storage_Žilinský | Žilinský | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 6.0 | 0.9600 | 0.96 | 0.0 | 0.0 | 0.0 | |
| hydrogen_storage_Banskobystrický | Banskobystrický | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 336.0 | 0.6532 | 0.50 | 0.0 | 0.0 | 0.0 | |
| hydrogen_storage_Bratislavský | Bratislavský | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 336.0 | 0.6532 | 0.50 | 0.0 | 0.0 | 0.0 | |
| hydrogen_storage_Košický | Košický | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 336.0 | 0.6532 | 0.50 | 0.0 | 0.0 | 0.0 | |
| hydrogen_storage_Nitriansky | Nitriansky | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 336.0 | 0.6532 | 0.50 | 0.0 | 0.0 | 0.0 | |
| hydrogen_storage_Prešovský | Prešovský | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 336.0 | 0.6532 | 0.50 | 0.0 | 0.0 | 0.0 | |
| hydrogen_storage_Trenčiansky | Trenčiansky | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 336.0 | 0.6532 | 0.50 | 0.0 | 0.0 | 0.0 | |
| hydrogen_storage_Trnavský | Trnavský | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 336.0 | 0.6532 | 0.50 | 0.0 | 0.0 | 0.0 | |
| hydrogen_storage_Žilinský | Žilinský | PQ | 0.0 | True | 0.0 | inf | -1.0 | 1.0 | 0.0 | ... | False | NaN | True | True | 336.0 | 0.6532 | 0.50 | 0.0 | 0.0 | 0.0 |
16 rows × 29 columns
In [28]:
# Adding shedding generator:
n.madd(
"Generator",
n.buses.index,
suffix=" load shedding",
bus=n.buses.index,
carrier="load",
sign=1e-3, # Adjust sign to measure p and p_nom in kW instead of MW
marginal_cost=1000, # Eur/kWh
p_nom=1e9, # kW
)
Out[28]:
Index(['Banskobystrický load shedding', 'Bratislavský load shedding',
'Košický load shedding', 'Nitriansky load shedding',
'Prešovský load shedding', 'Trenčiansky load shedding',
'Trnavský load shedding', 'Žilinský load shedding'],
dtype='object', name='Bus')
1) Solving The Model without a limit on CO2 emissions:¶
In [29]:
# Calculating the CO2 emissions for each region based on the resulted power output:
e = (
n.generators_t.p
/ n.generators.efficiency
* n.generators.carrier.map(n.carriers.co2_emissions)
)
e.tail(20)
# Calculating the CO2 emissions in total (for the whole country)[tonnes/MWh thermal]:
display(e.sum().sum())
# Adding global constraint for CO2 emission
n.add(
"GlobalConstraint",
"emission_limit",
carrier_attribute="co2_emissions",
sense="<=",
constant=e.sum().sum() * 1, # no reduction constraints in CO2 emission
)
# Solving the Model:
n.optimize(solver_name="highs")
n.generators_t.p # Retrieveing the power output
0.0
INFO:linopy.model: Solve problem using Highs solver INFO:linopy.io:Writing objective. Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:13<00:00, 1.08it/s] Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:02<00:00, 2.52it/s] INFO:linopy.io: Writing time: 17.16s INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log. INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 954856 primals, 2119937 duals Objective: 4.59e+08 Solver model: available Solver message: optimal INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
Out[29]:
| Generator | Coal_Košický | Coal_Trenčiansky | Gas_Bratislavský | Gas_Nitriansky | Gas_Trnavský | Hydro_Banskobystrický | Hydro_Bratislavský | Hydro_Košický | Hydro_Trenčiansky | Hydro_Trnavský | ... | Solar_Žilinský | Wind_Žilinský | Banskobystrický load shedding | Bratislavský load shedding | Košický load shedding | Nitriansky load shedding | Prešovský load shedding | Trenčiansky load shedding | Trnavský load shedding | Žilinský load shedding |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| snapshot | |||||||||||||||||||||
| 2017-01-01 00:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | 247.155251 | ... | -0.0 | 271.806237 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-01-01 01:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | 247.155251 | ... | -0.0 | 295.034377 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-01-01 02:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | 247.155251 | ... | -0.0 | 313.706895 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-01-01 03:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | 247.155251 | ... | -0.0 | 362.456009 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-01-01 04:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | 247.155251 | ... | -0.0 | 423.834183 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2017-12-31 19:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | -0.000000 | 18.437215 | -0.000000 | -0.000000 | ... | -0.0 | -0.000000 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-12-31 20:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | -0.000000 | ... | -0.0 | 485.519601 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-12-31 21:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | -0.000000 | ... | -0.0 | -0.000000 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-12-31 22:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | -0.000000 | ... | -0.0 | -0.000000 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-12-31 23:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | -0.000000 | 18.437215 | 6.369863 | -0.000000 | ... | -0.0 | 414.539877 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
8760 rows × 37 columns
In [30]:
# Regional distribution of loads in Slovakia
load = n.loads_t.p_set.sum(axis=0).groupby(n.loads.bus).sum()
display(load)
fig = plt.figure(figsize=(15, 6))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.EqualEarth())
# Plotting the network with bus sizes proportional to the load
n.plot(
ax=ax,
bus_sizes=load / 2e9,
geomap=True,
color_geomap=True,
projection=ccrs.EqualEarth()
)
svk_regions.plot(ax=ax, facecolor='lightblue', edgecolor='black', linewidth=1.5, transform=ccrs.PlateCarree())
ax.set_extent([minx, maxx, miny, maxy], crs=ccrs.PlateCarree())
plt.title('Regional Distribution of Load')
plt.show()
bus Banskobystrický 3.480420e+06 Bratislavský 4.103476e+06 Košický 4.391560e+06 Nitriansky 3.778554e+06 Prešovský 4.552601e+06 Trenčiansky 3.215057e+06 Trnavský 3.186314e+06 Žilinský 3.876638e+06 dtype: float64
In [31]:
carrier_to_exclude = 'load' # Excluding the "load" carriers, which are created during the Shedding Generator defining
load_excluded_generators = n.generators[n.generators.carrier != carrier_to_exclude]
dispatch_by_carrier = n.generators_t.p.groupby(load_excluded_generators.carrier, axis=1).sum().div(1e3)
#display(p_by_carrier.head(20))
fig, ax = plt.subplots(figsize=(15, 8))
dispatch_by_carrier.iloc[:720].plot(
kind="area",
ax=ax,
linewidth=0,
cmap="tab20b",
)
ax.legend(ncol=5, loc="upper left", frameon=False)
ax.set_ylabel("GW")
ax.set_title('Hourly dispatch per technology- JAN 2017')
ax.set_ylim(0, 9);
#########################################################
fig, ax = plt.subplots(figsize=(15, 8))
dispatch_by_carrier.iloc[:24].plot(
kind="area",
ax=ax,
linewidth=0,
cmap="tab20b",
)
ax.legend(ncol=5, loc="upper left", frameon=False)
ax.set_ylabel("GW")
ax.set_title('Hourly dispatch per technology through a single day- 1.JAN 2017')
ax.set_ylim(0, 9);
In [32]:
# The aggregate dispatch of the pumped hydro storage units and the state of charge throughout the day:
fig, ax = plt.subplots(figsize=(15, 8))
battery_storage_units = n.storage_units[n.storage_units.carrier == 'battery storage']
p_battery_storage = n.storage_units_t.p[battery_storage_units.index].sum(axis=1).div(1e3)
state_of_charge_battery = n.storage_units_t.state_of_charge[battery_storage_units.index].sum(axis=1).div(1e3)
hydro_storage_units = n.storage_units[n.storage_units.carrier == 'hydrogen storage underground']
p_hydro_storage = n.storage_units_t.p[hydro_storage_units.index].sum(axis=1).div(1e3)
state_of_charge_hydro = n.storage_units_t.state_of_charge[hydro_storage_units.index].sum(axis=1).div(1e3)
total_storage_region1 = p_hydro_storage + p_battery_storage
total_SOC_region1 = state_of_charge_hydro + state_of_charge_battery
#display(p_hydro_storage.iloc[:24])
total_storage_region1.iloc[:24].plot(label=" dispatch [GW]", ax=ax)
total_SOC_region1.iloc[:24].plot(label=" state of charge [GWh]", ax=ax)
ax.grid()
ax.legend()
ax.set_ylabel("MWh or MW")
ax.set_title('Total Storage: Dispatch and State of Charge throughout the day in Banskobystrický ')
plt.show()
###############################################################
# The aggregate dispatch of the pumped hydro storage units and the state of charge throughout a Month (January):
fig, ax = plt.subplots(figsize=(15, 8))
total_storage_region1.iloc[:720].plot(label=" dispatch [GW]", ax=ax)
total_SOC_region1.iloc[:720].plot(label=" state of charge [GWh]", ax=ax)
ax.grid()
ax.legend()
ax.set_ylabel("MWh or MW")
ax.set_title('Total Storage: Dispatch and State of Charge throughout a month in Banskobystrický ')
plt.show()
##############################################################
# the aggregate dispatch of the pumped hydro storage units and the state of charge throughout the whole year:
fig, ax = plt.subplots(figsize=(15, 8))
total_storage_region1.plot(label=" dispatch [GW]", ax=ax)
total_SOC_region1.plot(label=" state of charge [GWh]", ax=ax)
ax.grid()
ax.legend()
ax.set_ylabel("MWh or MW")
ax.set_title('Total Storage: Dispatch and State of Charge throughout the whole year in Banskobystrický ')
plt.show()
In [33]:
# Locational marginal prices:
fig = plt.figure(figsize=(7, 7))
ax = plt.axes(projection=ccrs.EqualEarth())
norm = plt.Normalize(vmin=0, vmax=50) # €/MWh
n.plot(
ax=ax,
bus_colors=n.buses_t.marginal_price.mean(),
bus_cmap="plasma",
bus_norm=norm,
)
cax = fig.add_axes([0.95, 0.15, 0.02, 0.7])
cb = plt.colorbar(
plt.cm.ScalarMappable(cmap="plasma", norm=norm),
label="LMP [€/MWh]",
shrink=0.6,
cax=cax,
)
ax.set_extent([minx, maxx, miny, maxy], crs=ccrs.PlateCarree())
plt.show()
In [34]:
display(n.generators['marginal_cost'])
generation_cost = (n.generators_t.p * n.generators['marginal_cost']).sum(axis=1)
sorted_costs = generation_cost[np.argsort(-generation_cost)] # Sorting generation costs in ascending order
cumulative_probability = np.arange(1, len(sorted_costs) + 1) / len(sorted_costs) # Calculating the cumulative probability
# Plotting:
plt.figure(figsize=(10, 6))
plt.plot(cumulative_probability * 100, sorted_costs, linestyle='-')
plt.xlabel('Cumulative Probability (%)')
plt.ylabel(' price (€)')
plt.title('Cost Duration Curve')
plt.grid(True)
plt.show()
Generator Coal_Košický 31.431739 Coal_Trenčiansky 31.431739 Gas_Bratislavský 44.520508 Gas_Nitriansky 44.520508 Gas_Trnavský 44.520508 Hydro_Banskobystrický 0.000000 Hydro_Bratislavský 0.000000 Hydro_Košický 0.000000 Hydro_Trenčiansky 0.000000 Hydro_Trnavský 0.000000 Hydro_Žilinský 0.000000 Nuclear_Nitriansky 13.778923 Nuclear_Trnavský 13.778923 Solar_Banskobystrický 0.010000 Wind_Banskobystrický 0.000000 Solar_Bratislavský 0.010000 Wind_Bratislavský 0.000000 Solar_Košický 0.010000 Wind_Košický 0.000000 Solar_Nitriansky 0.010000 Wind_Nitriansky 0.000000 Solar_Prešovský 0.010000 Wind_Prešovský 0.000000 Solar_Trenčiansky 0.010000 Wind_Trenčiansky 0.000000 Solar_Trnavský 0.010000 Wind_Trnavský 0.000000 Solar_Žilinský 0.010000 Wind_Žilinský 0.000000 Banskobystrický load shedding 1000.000000 Bratislavský load shedding 1000.000000 Košický load shedding 1000.000000 Nitriansky load shedding 1000.000000 Prešovský load shedding 1000.000000 Trenčiansky load shedding 1000.000000 Trnavský load shedding 1000.000000 Žilinský load shedding 1000.000000 Name: marginal_cost, dtype: float64
In [35]:
n.global_constraints.mu
Out[35]:
GlobalConstraint emission_limit -28014.414082 Name: mu, dtype: float64
2) Solving the model with a CO2 emission reduction of 100% (no emissions):¶
In [36]:
#Setting a constraint to solver:
n.global_constraints.loc["emission_limit", "constant"] = 0
In [37]:
n.optimize(solver_name="highs")
INFO:linopy.model: Solve problem using Highs solver INFO:linopy.io:Writing objective. Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:15<00:00, 1.06s/it] Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:02<00:00, 2.76it/s] INFO:linopy.io: Writing time: 18.9s INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log. INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 954856 primals, 2119937 duals Objective: 4.59e+08 Solver model: available Solver message: optimal INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
Out[37]:
('ok', 'optimal')
In [38]:
n.generators_t.p
Out[38]:
| Generator | Coal_Košický | Coal_Trenčiansky | Gas_Bratislavský | Gas_Nitriansky | Gas_Trnavský | Hydro_Banskobystrický | Hydro_Bratislavský | Hydro_Košický | Hydro_Trenčiansky | Hydro_Trnavský | ... | Solar_Žilinský | Wind_Žilinský | Banskobystrický load shedding | Bratislavský load shedding | Košický load shedding | Nitriansky load shedding | Prešovský load shedding | Trenčiansky load shedding | Trnavský load shedding | Žilinský load shedding |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| snapshot | |||||||||||||||||||||
| 2017-01-01 00:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | 247.155251 | ... | -0.0 | 271.806237 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-01-01 01:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | 247.155251 | ... | -0.0 | 295.034377 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-01-01 02:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | 247.155251 | ... | -0.0 | 313.706895 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-01-01 03:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | 247.155251 | ... | -0.0 | 362.456009 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-01-01 04:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | 247.155251 | ... | -0.0 | 423.834183 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2017-12-31 19:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | -0.000000 | 18.437215 | -0.000000 | -0.000000 | ... | -0.0 | -0.000000 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-12-31 20:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | -0.000000 | ... | -0.0 | 485.519601 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-12-31 21:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | -0.000000 | ... | -0.0 | -0.000000 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-12-31 22:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | 14.923516 | 18.437215 | 6.369863 | -0.000000 | ... | -0.0 | -0.000000 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
| 2017-12-31 23:00:00 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | 14.381279 | -0.000000 | 18.437215 | 6.369863 | -0.000000 | ... | -0.0 | 414.539877 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 | -0.0 |
8760 rows × 37 columns
In [39]:
carrier_to_exclude = 'load' # Excluding the "load" carriers, which are created during the Shedding Generator defining
load_excluded_generators = n.generators[n.generators.carrier != carrier_to_exclude]
dispatch_by_carrier = n.generators_t.p.groupby(load_excluded_generators.carrier, axis=1).sum().div(1e3)
#display(p_by_carrier.head(20))
fig, ax = plt.subplots(figsize=(15, 8))
dispatch_by_carrier.iloc[:720].plot(
kind="area",
ax=ax,
linewidth=0,
cmap="tab20b",
)
ax.legend(ncol=5, loc="upper left", frameon=False)
ax.set_ylabel("GW")
ax.set_title('Hourly dispatch per technology- JAN 2017')
ax.set_ylim(0, 9);
#########################################################
fig, ax = plt.subplots(figsize=(15, 8))
dispatch_by_carrier.iloc[:24].plot(
kind="area",
ax=ax,
linewidth=0,
cmap="tab20b",
)
ax.legend(ncol=5, loc="upper left", frameon=False)
ax.set_ylabel("GW")
ax.set_title('Hourly dispatch per technology through a single day- 1.JAN 2017')
ax.set_ylim(0, 9);
In [40]:
generation_cost = (n.generators_t.p * n.generators['marginal_cost']).sum(axis=1)
sorted_costs = generation_cost[np.argsort(-generation_cost)] # Sorting generation costs in ascending order
cumulative_probability = np.arange(1, len(sorted_costs) + 1) / len(sorted_costs) # Calculating the cumulative probability
# Ploting:
plt.figure(figsize=(10, 6))
plt.plot(cumulative_probability * 100, sorted_costs, linestyle='-')
plt.xlabel('Cumulative Probability (%)')
plt.ylabel(' price (€)')
plt.title('Cost Duration Curve')
plt.grid(True)
plt.show()
In [41]:
n.global_constraints.mu
Out[41]:
GlobalConstraint emission_limit -28014.414082 Name: mu, dtype: float64
In [42]:
# Locational marginal prices:
fig = plt.figure(figsize=(7, 7))
ax = plt.axes(projection=ccrs.EqualEarth())
norm = plt.Normalize(vmin=0, vmax=50) # €/MWh
n.plot(
ax=ax,
bus_colors=n.buses_t.marginal_price.mean(),
bus_cmap="plasma",
bus_norm=norm,
)
cax = fig.add_axes([0.95, 0.15, 0.02, 0.7])
cb = plt.colorbar(
plt.cm.ScalarMappable(cmap="plasma", norm=norm),
label="LMP [€/MWh]",
shrink=0.6,
cax=cax,
)
ax.set_extent([minx, maxx, miny, maxy], crs=ccrs.PlateCarree())
plt.show()
In [43]:
# Reduction wind and solar potential in 25% steps:
for reduction_step in np.arange(0.75, -0.01, -0.25):
previou_step = reduction_step + 0.25
reduction_factor = reduction_step / previou_step
for generator_name, generator_data in n.generators.iterrows():
if generator_data['carrier'] in ['Solar', 'Wind']:
n.generators.at[generator_name, 'p_nom'] *= reduction_factor
n.optimize(solver_name="highs")
print(f'{(1- round(reduction_step * 100))} % Reduction of Renewable Energy installable potential')
n.optimize(solver_name="highs")
print()
print(f'CO2 Shadow price in case of {(1- round(reduction_step * 100))} % Reduction of Renewable Energy installable potential:')
print()
display(n.global_constraints.mu)
print()
###################################################################################################
# Dispatch power per technology :
carrier_to_exclude = 'load' # Excluding the "load" carriers, which are created during the Shedding Generator defining
load_excluded_generators = n.generators[n.generators.carrier != carrier_to_exclude]
dispatch_by_carrier = n.generators_t.p.groupby(load_excluded_generators.carrier, axis=1).sum().div(1e3)
fig, ax = plt.subplots(figsize=(15, 8))
dispatch_by_carrier.iloc[:720].plot(
kind="area",
ax=ax,
linewidth=0,
cmap="tab20b",
)
ax.legend(ncol=5, loc="upper left", frameon=False)
ax.set_ylabel("GW")
ax.set_title(f'Hourly dispatch per technology ({(1- round(reduction_step * 100))} % Reduction of RE potential)')
ax.set_ylim(0, 9);
#-----------------------------------------------------------
fig, ax = plt.subplots(figsize=(15, 8))
dispatch_by_carrier.iloc[:24].plot(
kind="area",
ax=ax,
linewidth=0,
cmap="tab20b",
)
ax.legend(ncol=5, loc="upper left", frameon=False)
ax.set_ylabel("GW")
ax.set_title(f'Hourly dispatch per technology through a single day ({(1- round(reduction_step * 100))} % Reduction of RE potential)')
ax.set_ylim(0, 9);
#######################################################################################
generation_cost = (n.generators_t.p * n.generators['marginal_cost']).sum(axis=1)
sorted_costs = generation_cost[np.argsort(-generation_cost)]
cumulative_probability = np.arange(1, len(sorted_costs) + 1) / len(sorted_costs)
plt.figure(figsize=(10, 6))
plt.plot(cumulative_probability * 100, sorted_costs, linestyle='-')
plt.xlabel('Cumulative Probability (%)')
plt.ylabel(' price (€)')
plt.title(f'Cost Duration Curve ({(1- round(reduction_step * 100))} % Reduction of RE potential)')
plt.grid(True)
plt.show()
##############################################################################
fig = plt.figure(figsize=(7, 7))
ax = plt.axes(projection=ccrs.EqualEarth())
norm = plt.Normalize(vmin=0, vmax=50) # €/MWh
n.plot(
ax=ax,
bus_colors=n.buses_t.marginal_price.mean(),
bus_cmap="plasma",
bus_norm=norm,
)
cax = fig.add_axes([0.95, 0.15, 0.02, 0.7])
cb = plt.colorbar(
plt.cm.ScalarMappable(cmap="plasma", norm=norm),
label="LMP [€/MWh]",
shrink=0.6,
cax=cax,
)
ax.set_extent([minx, maxx, miny, maxy], crs=ccrs.PlateCarree())
ax.set_title(f'Electricity price per region ({(1- round(reduction_step * 100))} % Reduction of RE potential)')
plt.show()
INFO:linopy.model: Solve problem using Highs solver INFO:linopy.io:Writing objective. Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:13<00:00, 1.10it/s] Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:02<00:00, 2.43it/s] INFO:linopy.io: Writing time: 17.02s INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log. INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 954856 primals, 2119937 duals Objective: 4.84e+08 Solver model: available Solver message: optimal INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
-74 % Reduction of Renewable Energy installable potential
INFO:linopy.model: Solve problem using Highs solver INFO:linopy.io:Writing objective. Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:19<00:00, 1.27s/it] Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:03<00:00, 1.77it/s] INFO:linopy.io: Writing time: 23.68s INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log. INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 954856 primals, 2119937 duals Objective: 4.84e+08 Solver model: available Solver message: optimal INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
CO2 Shadow price in case of -74 % Reduction of Renewable Energy installable potential:
GlobalConstraint emission_limit -28016.106241 Name: mu, dtype: float64
INFO:linopy.model: Solve problem using Highs solver INFO:linopy.io:Writing objective. Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:41<00:00, 2.74s/it] Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:08<00:00, 1.17s/it] INFO:linopy.io: Writing time: 50.74s INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log. INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 954856 primals, 2119937 duals Objective: 5.32e+08 Solver model: available Solver message: optimal INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
-49 % Reduction of Renewable Energy installable potential
INFO:linopy.model: Solve problem using Highs solver INFO:linopy.io:Writing objective. Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:39<00:00, 2.64s/it] Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:07<00:00, 1.08s/it] INFO:linopy.io: Writing time: 48.54s INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log. INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 954856 primals, 2119937 duals Objective: 5.32e+08 Solver model: available Solver message: optimal INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
CO2 Shadow price in case of -49 % Reduction of Renewable Energy installable potential:
GlobalConstraint emission_limit -19881.836628 Name: mu, dtype: float64
INFO:linopy.model: Solve problem using Highs solver INFO:linopy.io:Writing objective. Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:39<00:00, 2.65s/it] Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:08<00:00, 1.15s/it] INFO:linopy.io: Writing time: 49.3s INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log. INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 954856 primals, 2119937 duals Objective: 9.20e+08 Solver model: available Solver message: optimal INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
-24 % Reduction of Renewable Energy installable potential
INFO:linopy.model: Solve problem using Highs solver INFO:linopy.io:Writing objective. Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:40<00:00, 2.69s/it] Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:07<00:00, 1.08s/it] INFO:linopy.io: Writing time: 49.45s INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log. INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 954856 primals, 2119937 duals Objective: 9.20e+08 Solver model: available Solver message: optimal INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
CO2 Shadow price in case of -24 % Reduction of Renewable Energy installable potential:
GlobalConstraint emission_limit -7650.545927 Name: mu, dtype: float64
INFO:linopy.model: Solve problem using Highs solver INFO:linopy.io:Writing objective. Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:39<00:00, 2.65s/it] Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:07<00:00, 1.11s/it] INFO:linopy.io: Writing time: 49.0s INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log. INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 954856 primals, 2119937 duals Objective: 2.09e+13 Solver model: available Solver message: optimal INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
1 % Reduction of Renewable Energy installable potential
INFO:linopy.model: Solve problem using Highs solver INFO:linopy.io:Writing objective. Writing constraints.: 100%|████████████████████████████████████████████████████████████| 15/15 [00:40<00:00, 2.70s/it] Writing continuous variables.: 100%|█████████████████████████████████████████████████████| 7/7 [00:07<00:00, 1.14s/it] INFO:linopy.io: Writing time: 49.95s INFO:linopy.solvers:Log file at C:\Users\Mahtab\AppData\Local\Temp\highs.log. INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 954856 primals, 2119937 duals Objective: 2.09e+13 Solver model: available Solver message: optimal INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
CO2 Shadow price in case of 1 % Reduction of Renewable Energy installable potential:
GlobalConstraint emission_limit -4.999777e+06 Name: mu, dtype: float64
In [ ]: